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Abstract 

We study a family of generalized slope limiters in two dimensions for Runge-Kutta discon- 
tinuous Galerkin (RKDG) solutions of advection-diffusion systems. We analyze the numerical 
behavior of these limiters applied to a pair of model problems, comparing the error of the 
approximate solutions, and discuss each limiter's advantages and disadvantages. We then in- 
troduce a series of coupled p-enrichment schemes that may be used as standalone dynamic 
p-enrichment strategies, or may be augmented via any in the family of variable-in-p slope lim- 
iters presented. 
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§1 Introduction 

Generally when solving advection-diffusion equations — which are not strictly diffusion dominated 
— by way of, for example, finite element or finite volume techniques, one observes the presence 
of spurious oscillations in the solution space often brought about by the existence of shocks in the 
space of approximate solutions as well as from the presence of sharp and/or discontinuous profiles 
in the physical domain itself. Such ill-behaved approximate solutions have led to the development 
of numerous methods designed with the intent to consistently stabilize and "limit" the solution in 
order to deal with these oscillations, as they are seen to arise quite frequently in common scientific 
applications. For example, slope limiters are known to be of central importance in storm surge 
modeling [HI [32] in order to obtain, for example, well-behaved solutions in the presence of compli- 
cated free-boundary conditions along adapting shorelines. Limiting regimes are also of substantial 
importance in quantum hydrodynamic systems [U [29] and surface wave models [26] where they 
are used to reduce the oscillations caused by mathematical dispersion terms {i.e. nonlinear third 
order spatial derivative terms) that pervade, for example, tunneling solutions. In fact, slope limiters 
are of fundamental importance in most applications in standard fiuid dynamics, being employed 
commonly in compressible Navier-Stokes [31], Eulers02], and magnetofiuid [TTl [M] applications, 
not to mention the important role limiters play in the study of radiative transfer [T^] and kinetic 
theory [18]; just to note a handful. 

From a numerical perspective, it is clear that one should desire that even shock dominated 
solutions, like both their smooth and non-limited counterparts, converge in p as p Pmax- However, 
such convergence is fundamentally coupled to the behavior of the error accumulation with respect 
to one's chosen slope limiting methodology, which, it turns out, must operate over a larger number 
of degrees of freedom, respectively, as p increases. For example, in a hierarchical basis (as shown 
explicitly below) the degrees of freedom grow nonlinearly as a function of p and each degree of 
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freedom ends up carrying information of potentially pathological (or undesirable) overshoots and 
undershoots which have developed over the native (or non-limited) solution space. It turns out that 
this complication introduces a substantial technical difficulty in practice, which many papers on 
numerical shock capturing [U El [71 [121 CH [IS [251 [23, [28] tend to avoid addressing directly. Most 
noteworthy is the observation that slopelimiters tend to limit the coefficients in their chosen basis 
independently of each other, in the sense that each component is adjusted based on information 
about the surrounding solution on a relatively local submanifold of the total domain. It then follows 
directly that the application of the limiter grows nonlinearly in each timestep as a function of p. 
Since a limiter de jure introduces error into the FEM solution space each time it operates on the 
FEM solution, more applications of it (iteratively) to the solution space should, as a general rule, 
lead to greater error accumulation assuming the first application always introduces approximately 
the same amount of error. In fact, this is what we observe in each of our limiters de facto. However, 
we offer an alternative approach to this problem below which is both highly efficient and consistent 
with the more general setting of /ip-adaptivity. 

It perhaps comes as no surprise that the same type of complications do not arise with respect to 
the mesh size h. That is, the convergence in /i as /i \ h^[^ tends to arise as a natural consequence 
of the usual h convergence, where convergence seems essentially guaranteed in most reasonable 
limiting regimes, while the order of convergence most certainly is not [27] . This issue raises another 
subtle technical difficulty which we will not address directly in this paper, though we will mention 
its importance in the proper context. 

Another important technicality pertaining to computational efficiency arises with respect to the 
well-known Courant-Friedrichs-Lewy (CFL) condition. In this setting the temporal discretization 
is (partially) bounded from above by the spatial discretization. That is, in order to reach a higher 
order accuracy at a fixed value of h one must project onto a higher order polynomial basis in p, thus 
reducing the admissible timestep At of the scheme — which obeys an inverse relation by virtue of 
the CFL condition: At oc 1/p as discussed in [38] . 

Since this p-dependence on the solution accuracy runs counter to the CFL restriction in a 
practical computational sense, substantial effort has been invested in developing "smart schemes" 
which in some way are able to "sense" the appropriate place {e.g. x & Q) within the solution 
domain to enrich the polynomial order p, while keeping other areas either unaffected or adaptively 
de-enriching areas of "less importance." The ultimate goal of these schemes is to attempt to 
substantially improve the computational efficiency of the numerical scheme without ceding notable 
accuracy in the solution. In fact, it is generally theoretically true that when one couples adaptive 
/i-refinement to p-enrichment {i.e. /ip-adaptivity) an exponential improvement in the convergence 
scaling of the solution may be obtained [lOj. However, dynamic adaptive /i-refinement is beyond 
the current scope of this paper and will be addressed elsewhere. 

On the other hand several different schemes have been developed for dynamic p-enrichment of so- 
lutions (independent of /i- refinement), though many suffer the added complexity of being extremely 
system (PDF) dependent. The advantage of system dependent regimes is that such schemes often 
display very close coupling to the physics of the solution {e.g. energy methods as discussed in [HO]). 
The disadvantage is, of course, that the scheme is very system dependent and hence whenever a 
variable is added or changed the entire scheme must be recalculated; which is particularly trouble- 
some for systems of equations which are not mathematically well-posed. Other schemes rely on — 
in the FEM setting for example — the generalized features of numerical variational solutions and 
as a consequence often depend strongly on a relatively large array of user defined constants. These 
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schemes are obviously quite attractive from a meta-application perspective, where being able to 
deal with generalizable systems displaying complicated initial-boundary data generates great allure 
in itself. In this paper we focus on the latter class of solutions, as we are interested in schemes 
which may apply to a large and generalized class of PDEs, without being bound, ab initio, to any 
one particular system of equations. 

Nevertheless, in the present paper we restrict ourselves to the class of discontinuous Galerkin 
finite element methods, where the underlying basis is chosen such as to signify a ubiety of discon- 
tinuous solutions - that is, we turn our focus in this paper to shock-dominated solutions. In this 
setting we are interested in the situation where continuously adaptive p-enrichment is coupled to an 
adapting-in-p slope limiting regime. We view this setting as very attractive, since the discontinuity 
sensors for p-adaptation schemes are well established [331 E] to be good sensors for slope limiting 
methodologies as well, where the p-enrichment leads to stability and efficiency of the scheme while 
the slope-limiting further stabilizes the presence of spurious oscillations emerging near pathological 
discontinuities as so approximated to order p. 

The outline of this paper is as follows. In §2 we present our generalized setting, which can be 
summarized as: given an advection-diffusion system of equations, consider the initial free boundary 
value problem recast into the weak formulation and spatially discretized. We then take a temporal 
discretization via a RKSSP DG approach in which we obtain the form of our approximate solutions. 
Our formulation is general, while our examples focus on problems of hyperbolic transport saving 
the more general applications for the sequel to this paper. In §3 we introduce a number of slope 
limiters consistent with any order p basis. The first is the vertex limiter regime of [25] , the second 
the classical Barth-Jespersen limiter [6j, and the third and fourth are minor adaptations of the 
former two limiters made with an eye towards improving the L^-error convergence by adjusting 
a so-called "blind spot" present in the previous schemes at higher order p. The fifth approach is 
comprised of a family of hierarchical reconstruction approaches [U [27], while the sixth regime is a 
linear restriction method that can be viewed as a generalization of a limiter originally sketched in 
[?]. The final limiting regime we present is a mixed extension of the previous limiters referred to 
here as a hierarchic recombination approach. Section §4 then provides numerical experimentation 
using the schemes presented in §3 - namely a classical advective scalar transport problem, and a 
stationary solution to a closely related problem with highly singular initial data. We also show some 
convergence results on an analytic test case. Finally, in §5 we present the adaptive p-enrichment 
schemes, which are fully coupled to the slope limiters from §4 ab initio. These come in two basic 
types, the first for (heuristically) smooth solutions, and the second for solutions demonstrating 
(vaguely) "appreciable gradients." 



§2 Advection-diffusion systems in the DG formalism 

We are interested in solutions to an initial-boundary value problem for a generalized advection- 
diffusion system of arbitrarily mixed hyperbolic-parabolic type in x (0,T), where f2 C with 
boundary dQ, such that the system satisfies: 

Ut + Fx — Gx = g, given initial conditions U\t=o = Uq, (2.1) 

and generalized componentwise Robin boundary values 



ttiUi + VxUi^x ipi ■ n + Ci ■ t) - fi = 0, on dVL. 



(2.2) 



^2 Advection-diffusion systems in the DG formalism 



That is, the system is comprised of a generahzed m-dimensional state vector U = U(t, x) = 
{Ui, . . . , Urn)-, an advective flux matrix F = F{U), a viscous flux matrix G = G{U, U^), and a 
source term g = g(t,x) = {gi, . . . , Qm), where as G and t G (0,T). The vectors a, b, c and 
/ are comprised of the m functions, = ai{t,x), hi = bi{t,x), Ci = Ci{t,x) and fi = fi{t,x) for 
i = 1, . . . ,m, where n denotes the unit outward pointing normal and r the unit tangent vector. 
In addition, because we are interested in approximate numerical solutions of the form of [2l |3] 



restricted in part to the family of methods for elliptic equations, we rewrite (2.1) as a coupled 
system in terms of an auxiliary variable S, such that 

Ut + F,-G, = g, and S = C/„ (2.3) 

where we have substituted in the viscous flux matrix the auxiliary term, so that G = G{U, S). 

For notational completeness we adopt the following discretization scheme motivated by [131 El] ■ 
Take an open i7 C with boundary dQ, given T > such that Qt = ((0,T) x Q). Let denote 
the partition of the closure of the polygonal triangulation of Q, which we denote Qh, into a finite 
number of polygonal elements denoted fie, such that ^ = {fiei, ^62? • • • ? ^en^}' ne E N the 
number of elements in Qh. In this work we define the mesh diameter h to satisfy h = minij{dij) for 
the distance function dij = d{xi, Xj) and elementwise edge vertices Xi, Xj G dQe when the mesh is 
structured and regular. For unstructured meshes we mean the average value of h over the mesh. 

Now, let Tij denote the edge shared by two neighboring elements Qe, and Qej, and for i G / C 
Z+ = {1, 2, . . .} define the indexing set r{i) = {j G / : flgj is a neighbor of fie,}- Let us denote all flei 
containing the boundary dflh by Sj and letting C Z~ = { — 1, —2, . . .} define s{i) = {j G Ib '■ Sj 
is an edge of flet} such that Tij = Sj for fle^ G flh when Sj G dfle,, j G Ib- Then for Sj = r{i) Us{i), 
we have 

dne,= [j r,j, and dne,ndnh= [j r^,-. 

We are interested in obtaining an approximate solution to U at time t on the finite dimensional 
space of discontinuous piecewise polynomial functions over fl restricted to given as 

= {v : v^n^^ G ^^(fi^J Vfi^, e ^4 

for ^^{flej the space of degree < p polynomials over fig.. 

Choosing a set of degree p polynomial basis functions Ni G ^^{flei) for i = 1, . . . ,np corre- 
sponding to the degree of freedom, we can denote the state vector at time t over fig,, by 

Up 

Uh{t,x) = Y,U\{t)Nl{x), ^xefle,, (2.4) 

£=1 

where the iV^*'s are the finite element shape functions in the DG setting, and the t/^'s correspond 
to the unknowns. We characterize the finite dimensional test functions 



Vt„u;t,eW'''^{Qh,^h), by Vf,{x) = J2'"'iNl{x) and u^^ix) = ^,^1^1 



X 



where v} and u:} are the coordinates of the test functions in each fie,, and with the broken Sobolev 
space over the partition ^ defined by 
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Thus, for U a classical solution to (2.3), multiplying by or uj^ and integrating elementwise 
by parts yields the coupled system: 



di 



/ U ■ Vhdx + {F ■ Vh)xdx - F : v^dx 
Jue. Jnc,- JUe 



S ■ ujhdx 



G ■ Vh)xdx + G : v^dx 
(U ■u;h)xdx- 



Vh ■ gdx, 



(2.5) 



U : iv'^dx = 0, 



where (:) denotes the scalar product. 

Now, let riij be the unit outward normal to dQsi on Tij, and let v\ri- and f |r ^ denote the values 
of V on Tij considered from the interior and the exterior of Qei, respectively. Then by choosing 
componentwise approximations in (2.5) by substituting in (2.4), we arrive with the approximate 
form of the first term of (2.5) given by, 

-f- [ Uh ■ Vhdx ^ -f- 
dt n dt 



U ■ Vfidx, 



(2.6) 



the second term using an inviscid numerical flux by 

^i{Uh\r,,,Uh\r,„Vh) = J2 / ^iUh\r,,,Uh\r,, 



,n,j) ■ Vh\T,,d'~ 



and the third term in (2.5) by. 



f ^ 



Fh : v'^dx ^ F : v'^dx. 



(2.7) 



(2.8) 



that 



Next we approximate the boundary viscous term of (|2.5|) using a generalized viscous flux ^ such 

r,,,^ii) ■ Vh\r,^dE 



j&ii)''^^^ 1=1 



(2.9) 



while the second viscous term is approximated by: 

=y^(S,, Vh) = [ Gh: v'^Jx ^ [ G : v'^Jx. 



Finally the source g term of (2.5) is given to satisfy 



Vh ■ gdx. 



(2.10) 



(2.11) 



^2 Advection-diffusion systems in the DG formalism 



For the auxihary equation in (2.5) we expand it such that the approximate solution satisfies: 



V / tj{Uh\r,^,Uh\r,,,uJh\r,^,nij)dE = 0, 



(2.12) 



j6H(j) »^ 



where, 



2 

EE / t/(C//.|r,,C7.|r,„a;.|r.„n.,)dS^ J] 5^ [ ^(C7) 



I ■ {nij)iUJh\r,jdE 



given U a generalized numerical flux, and where 



/ S/i ■ u^hdx ^ S ■ u^hdx, and / Uh ■ oj^dx ^ U : uj^dx. 

Combining the above approximations and setting, ^ = J2n^ e% denoting the inner 

product 

<■ bhdx, 



we arrive at our approximate solution to (2.3) as the pair of functions {Uh,^h} for all t G (0,T) 
satisfying: 

The Discontinuous Galerkin formulation 



a) UheC\iO,Ty,Sl), S,gS^, 

b) ^{Uh, Vhhg + ^(C/h, Vh) - &iUh, Vh) 



(2.13) 



- ^(S/,, Uh, Vh) + ^CEh, Uh, Vh) = J^ivh, Xh, t), 

c) ^iU,i:h,Uh,u;h,u;''J=0, 

d) Uh{Q) = n^c/o, 

where Yih is a projection operator onto the space of discontinuous piecewise polynomials S^, and 
where below we always utilize a standard L^-projection, given for a function /q G L^(f2eJ such 
that our approximate projection /q G L^(f2eJ is obtained by solving, L /g h^hdx = L f^Vhdx. 
We provide several explicit simplifled examples of this generalized formalism below, though in the 
followup paper we address models motivated by more complicated dynamics {e.g. see [8| [20] - [22ll2l] ) 



that employ the full system of (2.13) including multicomponent reaction-advection-diffusion and free 
boundary conditions, etc. 



The discretization in time follows now directly from (2.13), where we employ a family of SSP 



(strong stability preserving, or often "total variation diminishing (TVD)") Runge-Kutta schemes 
as discussed in [361 137]. That is, for the generalized SSP Runge-Kutta scheme we rewrite ( 2.13[ )) in 
the form: Mt/^ = R, where U = (C/i, . . . , Up) for each element from (2.4), where R = R(C/, S) 



8 



Adaptable RKDG systems 



is the advection-diffusion contribution along with the source term, and where M is the usual mass 
matrix. Then the generahzed s stage of order 7 SSP Runge-Kutta method (denoted SSP(s,7) or 
RKSSP(s,7)) may be written to satisfy: 

1/(0) ^ Jjn^ 

i-1 

C/(^) = J2 [oirrU' + AtA,M-^R") , for 2 = 1, . . . , s (2.14) 

r=0 

Jjn+l ^ jj(s)^ 

where R'' = R {U^, S'^, x, + 6r^t) and the solution at the n-th timestep is given as = U\t=t" 
and at the n-th plus first timestep by U"'^'^ = U\t=t"+^, with t''^'^^ = + At. The air and 
are the coefficients arising from the Butcher Tableau, and the third argument in R'" corresponds to 
the time-lag complication arising in the constraints of the TVD formalism. That is 6r = ^^1=0 f^ri, 
where /ij^ = Pir + Sz=r+i f^irC^ih where we have taken that air > satisfying X]r=o '^ir = 1- 



It is often possible to optimize the generalized SSP schemes of (2.14) by restricting to an op- 
timization class of stage exceeding order SSP Runge-Kutta time discretizations of [21] as long as 
p < 3. This class of SSP Runge-Kutta schemes has the advantage of optimizing the polynomial 
order p of the approximate solution Uh with respect to the r stage of the SSP Runge-Kutta scheme 
(incidentally satisfying SSP(r, p+ 1)) in order to minimize the effect of the rigid constraint intro- 
duced by the CFL condition on the timestep At. The limitation on p (i.e. requiring p < 3) is 
generally more restrictive than we encounter here, and thus, as will become apparent below, in 
the context of dynamic p-enhched slope limited solutions we are generally unable to exploit these 
optimization schemes directly. 



§3 A dynamic— in— ^ family of slope limiters 
§3.1 A transformation of basis 

Finite element approximate solutions are recovered with respect to any number of different finite 
element bases {e.g. Legendre polynomials, Lagrange polynomials, Labotto polynomials, Jacobi poly- 
nomials, Gegenbauer polynomials, Chebyshev polynomials, Bernstein polynomials, Gram-Schmidt 
polynomials, NURBS, T-splines, Wachspress functions, etc.). As a consequence of this, it is often 
advantageous to develop a strategy to transform into a specific basis in order to limit the solution, 
and then to transform back into the native bases to perform the remainder of the calculations. This 
occurs because some slope limiting regimes use fundamental properties of a certain choice of basis 
in order to develop a limiting strategy. We provide an explicit example of this procedure below, 
in the case of transforming between the Dubiner basis and the Taylor basis; or as we denote it 
below: by way of the invertible Dubiner-Taylor transform C. We also not here that for the sake 
of providing explicit calculations, we restrict below to triangular meshes, though the formalism can 
be easily extended to a more general framework. 



Take a solution vector U with approximate form Uh ~ C/ as given by (2.4), and project it onto 
the degree p Dubiner basis such that: 

Uhix,t)in^= J2 Uij{t)(j)i,{x), Mxe^e, (3.1) 

0<i+i<p 
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where the 4>ij{x) are the Dubiner basis functions for each degree of freedom in the solution vector. 

It is our aim to take this approximate solution Uh and limit it with respect to the k-th order 
Taylor basis via, for example, the vertex slope limiter of and the hierarchical reconstruction 
of [U [27], etc. Now, the Taylor basis in two dimensions is given to arbitrary differential order 
^ ^ + j) by the Taylor polynomial centered at c via: 



Uhix,y) = U 



h\c 



0<i+j<k 



dx^dy^ 



where Xc and yc are explicitly chosen as the values at the centroid c = (xc, yc 
Vie in the physical space VL — that is, each f2e. taking coordinates x E — 
2 + j > 1 in the sum denotes the differential order of the basis expansion (z.e 

1,3 en). 

Now, for cell averages satisfying U 
written by 



(3.2) 

of each finite element 
where it is clear that 
the indices satisfy 



Uh{x,y) = Uh\c + 



E 

0<i+j<k 



Qe\ ^ Jq Uhdx, the average of (3.2) may be simply 

(3.3) 



{x - XcYiy - ycV 



dx^dy^ 



such that subtracting (3.3) from (3.2) formally yields 




XcYiy-ycV {x-XcYiy 



iiji 




(3.4) 



Additional analysis (also see [28]) has shown empirically that the conditioning of the system 
in the Taylor basis (with respect to, for example, the invertibility of the Taylor mass matrix) is 
improved by rescaling with respect to the cell averages over the local bounds, given by ipAx = 
{xraa.x — Xrain) and ipAy = (^max — l/min) where ip = p foT p > 2, and '?/' = 2 for p < 2. It is useful to 
note here that in the master element representation these scalings are merely a pair of constants, 
while in the physical element representation they will in general vary elementwise. 

Then we are interested in implementing a locally renormalized Taylor basis prescribed with 
respect to the physical space Q given componentwise via the explicit formulation: 



¥^ij{x,y) 



[X — Xr 



{y - VcY 



where again cell averages are chosen to satisfy, 

1 



iy - yc) 

jlAy^ 



(3.5) 



X 



ilAx^ 



{y - yc) 

jlAy'J 



I a 



ilAx"- 



{y - yc)^ 

j\Ay3 



dxdy. 



Notice also that the constant terms of (3.4) vanish with respect to the barycenter c, which is 



just to say that the value of the centroid is by definition the cell average. Moreover, note that the 
renormalization vanishes for linear terms, since the average value is achieved at the centroid c (see 
for more examples at order p <2). 



Now we see that (3.4) satisfies in vector form that: 
Uh = Uh^oo + 



J2 '^^i 

0<i+j<k 



dx^dy^ 



Ax'Ay^ 



(3.6) 
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where we have denoted our effective Taylor basis ipij G ]R[f2], such that <^ij = in the poly- 

nomial ring ]R[f2] such that x E Q. By the polynomial ring ]R[fi] we simply mean the set of all 
polynomials with coefficients in M centered at a particular x E Q. The bracketed terms in (3.6) 
here represent our effective scaled coefficients, and from here forward the scaling parameters will 
generally be suppressed for notational simplicity. 

We will further make use of the fact that (3.6) may be viewed as the fc-jet over M^. That is, 
for r2 C and components of the approximate solution vector Uh the Taylor basis functions ipij 
comprise the abstract indeterminates of the fc-jet {J^Uh){^Pij) centered at c, in that by definition 

U,\n^^ := iJ^U,)iif,,), (3.7) 

such that our approximate solutions are elements of the abstract jet space Uh\ne- ^ JciM-^, The 
jet space J^(]R^, Q) is simply defined as the set of equivalence classes of fc-jets which agree to order 
k and map between the Cartesian plane and an element of Q, as clearly our approximate solutions 
in the Taylor (polynomial) basis do. By the set of equivalence classes of fc-jets which agree to order 
k, we mean for any two solutions Vh\n^. and Uh\ne- in the Taylor basis restricted to flej — that is 
fc-jets — the equivalence relation Uh\n^. — Vh\ne. ~ holds to order k. 

In this sense, an effective slope limiter may be viewed as a stabilization rescaling of the jet by the 
k coefficients (as derived in §4), such that the slope limited approximate solution vector UJ^ is 

formally the same as the stabilized fc-jet centered at c; that is C/J^Iq^. := {J^oXJ h){Vij) where both 
the approximate solution and the corresponding limited approximate solution are each, respectively, 
elements of the same abstract jet space U h\Q.^ . ,UWq_^ , G J^{M?,VL) when letting the equivalence 
relation ~ be approximate ~/i {i.e. approximate with respect to the solution order accuracy but 
with vanishing asymptotics). 

Now, in order to work between the Taylor basis representation ^pij and the Dubiner basis rep- 
resentation we must construct a transformation between the physical element space VL and 
the master element space A^, as well as a transformation between the two (abstract) polynomial 
bases. Below we make these mappings explicit, and refer to them collectively in this work as the 
Dubiner-Taylor transform, which is given by the invertible mapping £: ]R[A^] — )■ J^{E?,VL). 

First consider the usual Dubiner basis functions in the master element space componentwise 
0jj G ]R[A^] for ipij = (j)ij{x), and M[A^] the polynomial ring in coordinates x E Ai given by: 

'P^,=prii^i)[^Jpr'M^2), (3.8) 

using p-th order Jacobi polynomials with weights a, (3, such that Pp'^{-) is evaluated with respect 
to the coordinates a; = (^, r^) of the master triangle element, where the master element quadrilateral 

transformation in the Dubiner mapping provides that: ipi = ^ "^{i-^) ~ "^2 = such that 

ipi = ipiix) and ip2 = 'ip2{x). 

Now, consider the two state vectors, = (0oo; 0iO) • • • > 0cd)"^ and (p = {ipQQ,(piQ, . . . ,ipcd)'^ , 
where in the lexicographic ordering (described in detail in §3.2) we have c + d < p. Now, we may 
transform between the master and physical element representations of our components (p = ip{x,y) 
and = (f){C,, T]) using the following affine mapping: 



X 



-^{^{xi-x2)+ vi^i - xs) - X2 - xa^ y = -^\^^{yi - y2)+v{yi - ys) - y2 - 



(3.9) 



3 A dynamic-in-p family of slope limiters 
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Figure 1: We look at the maps N: R[M] Jl:{R^,M), S: J^{R'^,^) J^^(M^A^), and 
C: M.[Ai] —7- J^(]R^,i7), where and J^m are the abstract operators that limit in either the 
physical element space Q or the master element space Ai. 



with inverse given by 



C = x{ (2/3 - yi)ix - ^{x2 + xs)) + {xi - X3)iy - ^{y2 + 2/3)) }, 

1 1 

V = X-{ {yi - y2){x - -{x2 + 0:3)) + {x2 - xi){y - -{y2 + 2/3)) 



(3.10) 



Here {(xi,|/i), {x2,y2), (2^3,2/3)} are the vertices of the triangles in the physical space, and the area 
of the physical element fie is given from the cross product of two of the triangle edge vectors, 
via the usual formula 



X = 2 (X22/3 - 2:32/2 + 3:32/1 - Xiy3 + Xiy2 - X22/1) 



-1 



Then by substitution of (3.9) and (3.10), we may easily construct the invertible mapping 
S: J^(]R^,fi) — )■ J^(]R^,A^), such that <j = S{(p) represents the Taylor basis in the master ele- 
ment space Ai. That is, to construct S explicitly we take the constant first order transformation 
rules for the derivatives in the base coordinates, given by 

dx^ = Xiy-i- yi), dy^ = xixi- X3), <9^?7 = x(2/i - 2/2), dyT] = x{x2- xi), (3.11) 

in the master element representation Clf,. G A4, and 

d^x = {x2 - Xi)/2, % = (1/2 - 2/i)/2, d^x = {xi - X3)/2, dr^y = [yi - y3)/2 (3.12) 

in the physical element representation Q^i ^ ^• 

Thus provided the coordinate pair rj) in the master element representation Qg. G we may 



use (3.9) evaluated at the element quadrature points i to fully determine S, where the evaluation 



at the quadrature points allows for explicit computation of the integral averages in the Taylor basis 



components (3.5), or, more explicitly, where we compute: 



(x - Xcy\ f (y - ycV 



i\Ax' 



jlAy^ 



' E^^^^^^^^^N^^^:^lldetJ| 



ilAxl 



jlAyj 
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for wi the quadrature weights and the determinant of the Jacobian matrix J satisfying |detJ| = 

I ftr 9j/ dx dy I 

All that remains then is to find the coefficient matrix which constructs the change of polynomial 
basis mapping N: ]R[A^] — )■ J^(]R^, Ai), such that we may write the components of the transformed 
Taylor basis <^jj, given by terms Tij^ij, with respect to the components of the master element frame 
Dubiner basis given by terms Dj^di^y, or such that we recover the matrices 



bij, given by terms Dy^.j, 

T = N(0), and vice versa D = N^"'^(<j). 



(3.13) 



But in light of (3.1) and (3.4) it follows that for the K-th component of the m-th size solution 
vector Uh in we may solve for the Taylor coefficients Tij using the system: 



\k^^''cdUd7]dij 



I In '^ood-vd^ In '^oo'^iodrid^ ... I^ <^oo'^cdd'nd^\ 

Id <;oo<iwdr]d^ 4 <;fodr]d^ ... ^w^cddrjdS, 

Gj 

yid^^c^QQ^cddTldi J^^^^ioq,ddrid^ ... Id^^c^l^diqdi j 



/Too\ 

Til 

\T.d) 



(3.14) 



for the K-th component of Uh- Note that for the convenience of the reader, we suppress the com- 
ponent index k here and below, though it should be understood that the slope limiting operations 
are generally performed componentwise over the elements of the solution state vector. 

Now, extending (3.14) over all the components, the Taylor mass matrix tensor on the right 
and the inner product matrix on the left serve to define the desired transformation: 

N(0) = M"^ oP^. 

Its inverse is simply given by forming the Dubiner mass matrix tensor M,^ and the inner product 
matrix in denoted P<^, such that: 

N(<^)-i = M^ioP^. 

Then we have fully constructed the invertible Dubiner-Taylor transform C: J^{E?^VL) as 

satisfying 

£(0) = S"^ o N = o o P^ = T o c^. (3.15) 

with inverse satisfying : 



C-\ip) = N(<^)-i o S(c^) = M , 1 o P^ o S(c^) = D o 0. 



§3.2 The formal vertex based hierarchical limiters 

We now formally construct the generalized vertex-based slope limiter based off the Barth-Jespersen 
limiter ^ |25] . In this context we define a neighborhood as comprised of those elements that share 
a common vertex Xi, indexed with respect to every vertex of each finite element cell VL^y More 
clearly, we define the focal neighborhood = {^ej}i (in the sense of the foci of geometric optics, 
as shown in Figure [2| as the collection of elements such that Xi G — where {l^e^ji includes the 
base element Qg. — such that i = 1, 2, 3 over triangular elements. 



3 A dynamic-in-p family of slope limiters 



13 




Figure 2: Here we show the focal neighborhood i7j of a base element filled in black. Green, 
red and blue are the three focal neighborhood groups based at vertices Xi of the black base cell, 
while purple are cells contained in more than one of the two focal neighbor stencils (incidentally 
comprising the edge neighborhoood of fiej- In a contrasting geometric locale, the orange base cell's 
edge neighbors are each filled in yellow, comprising the edge neighborhood VIe- See Figure |4] for 
details. 



We now note that one must choose a base space in which to implement this slope limiter {e.g. 
the physical Q or master Ai element spaces, etc.). A fairly common choice {viz. [251 |2B]) is to limit 
with respect to the full physical space Q. However, in the context of the local DG formulation this 
choice is not always so clearly taken. That is, given our transformations from §3.1, it is clear that 
we may not require the full Dubiner-Taylor transform C but rather have the option to restrict to 
the master element space Ai by simply using the invertible map N. More clearly, since local DG 
formulations often exploit computational efficiency by working over a master element representation 
Ai, we are presented with a choice of composition maps to limit in the master or physical element 
spaces as shown in Figure in and given either by N ^ o o N over , or by £ ^ o o £ over il. 
However, since (3.15) shows that C requires the extra algorithmic step of transforming back into the 
physical coordinate frame Q, in the name of computational efficiency, we clearly prefer the former 
composition given the context of a relatively standard local DG method. However, when working 
in a global DG formulation where one elects, for example, a global linear solve, it may be more 
beneficial to limit with respect to Q, which as shown in Figure [l] may also be easily accomplished. 

Now, we may define the explicit role of the vertex slope limiter as: a method of finding the limiter 
matrix ex. = (cki, . . . , CKm)"^ such that for the solution vector satisfying Uh = {Ui, . . . , Um)'^, with 
m the number of unknowns in the system of equations, a vector defined by ck = {a^^\ . . . , a'-'^^)^ 
for each order derivative i + j < k, the limiter coefficients a^'^'^^^ G [0, 1] allow for a recasting of 
the renormalized solution in (3.6) componentwise in the vertex slope limited form with respect to 
a focal stencil, that is C ^If for a fixed vertex Xi (see Figure |4] for more detail). 

In fact, regardless of the initial location containing the state vector {i.e. with respect to Ai or 
with respect to Q) by simply using our transformations S and N from §2 we can recast (3.6) in the 
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master element space Ai such that componentwise we have the vertex slope limited approximate 
solution which satisfies: 



0<i+j<k 



(3.16) 



where U and U correspond to the approximate solution vector Uh transformed to the master element 
frame in the Taylor basis representation. 

Now, notice that above there exists only one a*^*"'"-'^ for each top k-th order mixed derivative in 
^ and 7]. In order to recover the a*^*+-^^'s in the polynomial basis expansion, we must decompose our 
solution Taylor expansion into mixed order linear reconstructions. To do this, we first order our 
Taylor polynomial into a hierarchical basis such that each monomial index b = b{i,j) is provided 
using the lexicographic ordering with ordered lattice pairs given by the sequence (0,0) < 

(0, 1) < (1,0) < (0,2) < (1, 1) < . . . = (z,j) corresponding to indices b, respectively; that is by the 
sequence (1) < (2) < (3) < (4) < (5) < . . . < (6) . . . < (s) in the Taylor expansion. In fact, the 
monomial index in the hierarchy may be determined by the diophantine equation: 



|(j + l) + zj + ^(^ + 3) + l. 



(3.17) 



Then we generate the hierarchical triangular sequence s = s{p), where p = p{i-,j) satisfies p = (i+j), 
such that s determines the upper bound on the degrees of freedom in the polynomial expansion. 



1 



(p + l)(p + 2), given inverse g = g{s) such that g 



1 



(3.18) 



where [-J is the usual fioor function. Note that in particular we may use g = g{s) or g = g{b) for 



g{b) G i corresponding to level i ^ kop (defined below) since by virtue of the mapping (3.18) both 
return the same value. 

Then letting U^^^ ^ be the value of the 6-th term in the polynomial basis of Uh at the centroid c of 

element Vte^ containing Xi = {C,i,'rii) in the master element representation, we define the maximum 
^max minimum ?7™" values for each unknown monomial at Xi over the focal stencil Q /. situated 

with respect to the master element frame Qf. as 



rrmax 



max {U-^h,c} and Ul 



mm 
b 



mm 



(3.19) 



Now, we are able to define the {i + j)-th linear reconstructions U^^^''^ over the vertices Xi of 



any element by taking derivations with respect to the monomial coefficients of (3.16). That is, the 
linear perturbation of the constant term is constructed such that. 



(1) 



u 



dU,. 



da, 



drj 



iVi - Vc), 



for s 



(3.20) 



Moreover, it is now direct to construct the higher order terms whereby setting 



for b{i,j) > 1, 



3 A dynamic-in-p family of slope limiters 



15 



A schematic of the vertex-based methods 



i - the first index of the Taylor expansion 
j - the second index of the Taylor expansion 
b{i,j) - the monomial index 

- the Taylor monomial of index b 
^ht^^ - the linear reconstruction at {C,i,rji) 
U^^^, t/j™ - the extrema over the stencil flxi 
a^'^^ - the limiting coefficients on f2e, n 

142, m) 




Base element, Q^.^ G fi^. 



Centroid of Qe, G Qx 



(^3, ^3) 



Figure 3: Here we provide a key for the vertex and Barth-Jespersen limiters of §3.2-§3.2.1. Generally 
theses limiting procedures depend on the chosen stencil and a local linear reconstruction of the 



solution in order to develop the limiting coefficients from (3.16). 



such that for any mixed derivative order in the hierarchical basis — as a property of the lexicographic 
ordering — we can write: 



(i+J) 
bA 



(3.21) 



for any polynomial order k. Proceeding, we can now define the correction factors for each 

element figp where the vertex-based condition is simply defined as 



a 



mm < 



mm 



mm 



rrmax TT I 

0,1 i,c,b 



rrmin TT^l 
'^i,b '~'i,c,b 



i.c.b 



for a 



for a 



for a 



b,i 

b,i 



W 



i,c,l 



(3.22) 



which, again, is determined separately for each monomial represented in the master frame. 

These a^*^"''' determine a set of limiting constraints for every hierarchical monomial in the Taylor 
expansion, but as in |25], we minimize over derivatives of similar top order, such that we recover 
the components: 



a 



Kp) 



mm . 

9{b)=KPo) 



(3.23) 



Notice that these limiting coefficients span the level 1{pq), not the hierarchical index b corresponding 
to level l{p) (where p and po are fully explained below). That is, in the hierarchical basis the linear 
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Figure 4: Here we show the focal stencil Qf^ and the edge stencil fie. of a base element f2e, = A. 
The stencils are defined with respect to the base elements vertices Xi for i = 1,2, 3, such that the 
focal stencil at Xi is = {J, I, N, H, B, A}, and likewise Qf^ = {J, K, L, M,C, A} and Qf,^ = 
{C,G, F, E, D, B, A}. Similarly the edge stencils are given by: Qei = {JjBjA}, = {'J:C,A} 
and — {B,C,A}. Notice that the union of sets recovers the focal neighborhood {^f = ^i^fj 
and the edge neighborhood {flE = Uif^eJ, while the restriction of the symmetric difference of sets 
defines the focal neighborhood group (GjQ/. ) and edge neighborhood group (Gi^E. ) for any 
vertex j. 



reconstructions from the perturbation at the level below {i.e. level {I — 1)) are what effectively 
determine the limiting coefficient at level I {e.g. the gradient terms). More precisely, the level 
I — l(po) is determined with respect to the sequence of integers starting at po(po — l)/2 + 1 and 
increasing by one until reaching ^0(^*0 + l)/2- Then the level is defined by [ = sup{g{pQ{pQ — 
l)/2 + 1), . . . ,g{po{Po + l)/2)}, where Po = 1 for the strictly linear case, and in general is a positive 
integer such that po < p and is fully determined by g{b{i,j)). In general however, the level i{p) 
spans I = S'u.p{g{p{p + l)/2) + 1, . . . ,g{{p + l){p + 2)/2)} such that the level below ((po) simply 
corresponds to setting p — po — 1. 

Finally we limit the magnitude of the correction by the maximum value of every correction 
factor of greater than or equal order. In other words, we do not allow a higher order correction to 
demonstrate greater regularity than a lower order correction, and in fact empirical experimentation 
has found this to be a necessary constraint. That is, setting q = {i + j) and r — {i' + j') for i' and 
/ indices, then we determine an upper bound on the correction parameter by resetting: 

a^"^ := max \/q > 1, Vr > q. (3.24) 

q<r,l<ltop 

The top level (top simply corresponds to the level whose upper bound is determined by g{s) — g{b). 
Also notice that the derivative order {i + j) is fundamentally coupled to the level i, and so is in 
some ways redundant notation which we have used in order to emphasize this coupling. 
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It is also worth noting, that as a consequence of the above construction we are now easily able 
to implement an arbitrarily higher-order extension of the Barth-Jespersen limiter [HIISS], where we 



may perform the exact steps as above, but simply exchange (3.19) with 



[/fr=.max {U:^ and t/^''^ = . min {[/^^J, (3.25) 



where ClEi is the edge stencil of at Xi in the master element representation — or the corre- 
sponding set of those physical elements sharing an edge with Q^j at vertex Xi such that the base 
element fig ^ (see Figure [2|). 



A schematic is provided in Figure [3] which is meant to simplify the notation and unify the 
basic principles underlying both the vertex and Barth-Jespersen limiters (as well as the adapted 
vertex-based limiters of §3.2.1). 



§3. 2a On adapted vertex based limiters 

Both the vertex limiter and the Barth-Jespersen limiter from §3.2 demonstrate a similar — though 



often times non-ideal — behavior. That is, notice that in both the definition of (3.19 ) and 3.25 ) that 



we have found a maximum or minimum with respect to a local neighborhood of the mesh. Hence, 



in either case, when we compute the limiting coefficients in (3.22) a local bound {e.g. (3.25)) is 



always achieved, even in the degenerate case of when U^^™ = U^^^. 



As a consequence of this, the numerator in the quotients of (3.22) vanish on elements admitting 
a local extremum, leading to persistent and excessive diffusivity {i.e. limiting a = at each such 
timestep) arising at all orders in each local extrema of the mesh, even when those extrema are 
neither spurious nor potentially unstable; and moreover, this behavior compounds in p since as 
p increases the number of degrees of freedom {i.e. monomials) in the solution which have local 
extrema also increases nonlinearly. 

This behavior over values of local extrema can become quite dominant depending on the mesh 
geometry. In particular, since the vertex-based limiter has a larger local neighborhood {i.e. the focal 
neighborhood) than the Barth-Jespersen limiter, in principle it should provide more information 
from which to glean a more accurate approximate local reconstruction. However, due to this 
"diffusivity," the larger local neighborhood actually lends itself towards increasing the nonlocality 
of the diffusive effects of the neighborhood-wise extrema as p Pmax, and hence in practice can 
actually precipitate greater diffusion in the vertex limiter than the native Barth-Jespersen limiter 
as p increases (up to the mesh geometry) . 



In order to reduce this so-called "blind diffusion" in both limiters we introduce a simple func- 
tional which attempts to treat a portion of this special case separately. That is, we simply replace 
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K22^ with: 



a 



(i+j) 



mm < 



mm 



mm 



mm 



mm 




for 

for 
for 
for 



Tjii+j) 



Tjei 



> 



Tjei 



u, 



(i+j) 



(3.26) 



, for U, 



ii+j) 



Uf 



where /max, /min G (0, 1) are constants used to hmit the rate at which the extrema diffuse (that is, 
reduce the rate at which error is introduced into the solution), and when /^ax = /min we denote 
them by fd- 

We find when setting fd = lwe generally get a very moderate improvement in the limiting error 



behavior of both the vertex and Barth-Jespersen limiters. Nevertheless, clearly (3.26) has only 



accounted partially for the degenerate local extrema cases, in particular it still fails to properly 
account for the case of = UJ^^""^, and the absolute value is used to account for the fact that 
the signs have not been separately controlled. We have developed strategies for adopting fixes for 
these issues into the limiter, but in general find even the augmented regimes to still demonstrate 
substantially more diffuse behavior than the restricted regime presented in §3.4, and so will suppress 
any further comment on the subject at present, simply noting that it is possible to improve upon 
the basic behavior of the limiter in p by developing selection strategies to deal with the many special 
cases which arise over solutions locally, and where alternatively one is often also able to improve 
the error behavior by tuning /max and /min- 

It should be additionally noted here that in [25] a mass lumping strategy is implemented with 
respect to the triangular meshes in order to prevent the formation of undershoots and overshoots 
caused in the presence of the non-orthogonal Taylor mass matrix. It was also demonstrated in |25j 
that this strategy can have a measurably beneficial effect on the error behavior for p < 2, and is 
thus clearly worth further examination at higher p. We will return to this issue briefly in §4 as a 
note of comparison between the implement at ional strategies. 



\3.3 The hierarchical reconstruction via MUSCL or ENO 



We now consider the hierarchical reconstruction scheme presented in P and [27]. Formally in this 



setting we simply take derivatives of (3.2) in the master element frame, and work locally over the 



averages and differences of these differential reconstructions. The method is presented as a two step 
process, where we start in step 1 at the highest order derivatives and work down to the lowest, with 
the caveat that the linear and constant terms are dealt with separately in step 2. 

Step 1. Starting at the top order coefficient k, a linearization of the {k — l)-st derivative of 
(3.2) is given by (3.21 ) in the Taylor basis foT i + j = k — 1. Here, however, we recover the entire 



higher order component including the nonlinear terms, so that we must employ our monomial index 



function b{i,j) given in (3.17). 
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That is, beginning at the top level l{k) ioi i + j = k we define the linear part as satisfying: 



^linear 



\fb e lik) A Vfie, e fix, 



(3.27) 



where X here and below may he f or E {i.e. the focal and the edge neighborhoods, respectively, 
as shown in both Figure |2] and Figure |4]), and where here and below A is the logical conjunction 
operator and V is the corresponding logical disjunction operator. 

At the lower (nonlinear) levels {i.e. the levels i such that 1(1) < I < i{k)) by expansion - after 
recovering the i and j indices of the base 6-th component — then setting b = b{i + i',j + j') and 
integrating locally over each cell in the neighborhood, we have that: 



+ [ I 5^ T7r7!^5(^-^c)*'(^-^c)^'|rf^C V6 < s A Vfie, G fix. (3.28) 



Likewise for each level I we integrate the higher order perturbative terms such that: 



{i+i) 



Vcy'{^ - ^cV >drid^, V6 < s A Vfie, e fix. (3.29) 



It is then these two averages which serve to limit the lev el i co mpon ents o f the Taylor basis by way 

(3.30) 



of the linear type average 

Ub^^\n, of the difference of (|3.28D with (|3.29D in each h 



^linear 



Now the linear terms of (3.30) will be used to determine the candidates for the updated values of 
the base cell flbase] which is to say that the {k — l)-st component of the k-th order jet {J'^oU h){'iij) 
is limited by filtering a set of candidates through a family of minmod functions, such that: 



Tjei 

^base 



minmod^^^^gf,^, (piZ'Ln^, 



(3.31) 



Notice that we may also choose to find candidates over restricted subsets of the full neighborhood 
fix in order to try and more effectively localize our limiting. For example, we may choose to find 
the minmod function over the local stencil fix^ centered about a vertex of the cell and then perform 
a different selection rule over that set of candidates; or, alternatively, we may compute the integral 



averages over the local stencil fix^ in (3.27)-(3.29) and then perform a minmod with respect to the 



full neighborhood fix. We have implemented and tested a number of these different regimes, and 
consider each of them in this paper to live under the general heading of "hierarchical reconstruction 



schemes," though for the sake of brevity we focus only on (3.31) below 



Note that we perform Step 1 for each level {{j) where j < k, and recursing down to the level 
corresponding to the [ associated to the quadratic components at p = 2; where first we limit the 



difference (3.30) across the neighborhood of a base element in order to reconstruct the values on the 



base cell proper. For these purposes, we employ the following set of minmod^ = $^ functions. The 
MUSCL reconstruction method relies on the function: 



mim U, 



* '-'linear 1^ 



maxj ( Ufj^_^^^^^Q^ j , if f/^ 



0, 



if Uj,"^^^ o > Vfie, e fix, 

^linear i^ ^e; 



O < Vfie, G fix, 

linear ' 

otherwise, 
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A schematic of the hierarchical reconstruction method 

• i - the first index of the Taylor expansion 

• j - the second index of the Taylor expansion 

• b,b - the monomial/hierarchical indices 

• - the Taylor monomial of index b 

• U^^^^ o - the candidates of the reconstruction 



• f^;, n ^ the higher order terms 




(6,^1) (^3,^3) 



Figure 5: Here we provide a key for the hierarchical reconstruction method developed in §3.3. The 
limiting procedure depends on the entire neighborhood the fully integrated solution, and a 
choice of minmod functions in order to reconstruct the limited form of the monomial coefficients on 
the base cell. 



while the ENO reconstruction is given by 



* '-'linear 1^ 



^linear 'i^Sj^ 



if U, 



^linear 1^ 



= mm 



(i+i) 



^linear i^e 



Additionally, following [27j, the minmod^ function may be set as a center bias scheme given by 

= $tn I n.^A.^'^f \ jj{i+3) \ (3.32) 

X I \ / ^ y ^linear -S^ei J ' j. / j ^linear i^e^ I ' ^ ' 

\ k=l / 

or the weighted ENO scheme. 



$1^ 




k=l 



(i+j) 

^linear t^ej^ 



(3.33) 



where in either case r is the total number of neighboring cells of the base cell ^e^ase^ e is a user 



defined constant, and (■) in both (3.32) and (3.33) is merely standard multiplication. It is known 



that setting e large helps to achieve the expected order of accuracy over triangular meshes. 

Step 2. Now we address the case of how to limit the solution with respect to the linear i + j = 1 
and constant i = j = cases. For the linear case, we simply choose to limit with respect to a subset 
of limiting regimes, including those in §3.2 §3.3 §3.4 and §3.5. We choose this, in particular, in order 
to electively replace the MUSCL and ENO schemes from Step 1, which are relatively speaking more 
diffuse in our experiments at level 1(1) than some other possible alternatives. 
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Finally, the constant terms at level i{0) are simply set equal to the average value on their base 
cell, U^oo\ni,ase order to enforce invariance of the cell averages. In other words, the constant terms 
remain unchanged. 

We should also note that recent improvements have been made in the context of hierarchical 
reconstructuion techniques. In particular, recent work has been done to extend the formalism above 
to include WENO-type linear reconstructions. That is, the WENO-type formalism of [35l HI] has 
been extended to the context of hierarchical reconstruction based limiter regimes in |13] specifically 
in order to address the fact that the MUSCL and ENO type approaches have been shown to often fail 
to give the desired order of accuracy on triangular meshes. These techniques rely on the conditioning 
of a local (fix-restricted ) recontruction matrix, and are beyond the scope of the present paper. We 
direct the interested reader to [13] . 



§3.4 On a dynamically adaptive linear restriction 

In this subsection we generalize and update a version of the BDS limiter that had its foundations 
initially seeded in [7J for linear polynomials over uniform structured meshes. We present the formal 
construction of a substantially more general form of this limiter, to act over an arbitrary order p 
basis by way of a linear restriction technique over unstructured triangular meshes. This limiter is 
developed with an eye towards p-enrichment schemes, and in particular /ip-adaptive schemes, where 
in areas of high (jump) variability one generally wants to reduce the order of p while refining the 
mesh parameter h. In this section we first restrict back to the Dubiner basis (pij G ]R[A^], in part 
to compare to the same implementation carried out in the Taylor basis as a consequence of the 
formulation presented in §3.5, which for linears turns out to be equivalent. 

Let us first restrict to the sub-quadratic terms of the basis for any order p, such that we are only 



concerned initially with the terms corresponding to i + j < 1. Then, similar to (3.19), setting U^^ 
as the constant piece of the Dubiner basis in Uh of the base element Clej containing Xi = {C,i,f]i) 
in the master element representation, we define the maximum JJf^'^^ and minimum [/P™ values for 
each unknown at every Xi G Clej over the chosen stencil Clx^ as 

f/f^"= max {U-'} and Uf''= min {U-'}. (3.34) 

Next we take the full approximate solution restricted to its sub-quadratic part and evaluated at 
the three vertices of the cell, denoted by the three values U{x^)\iJ^j<i for i = 1,2,3 corresponding 
to the vertices, while i + j corresponds to the polynomial order. Then at each vertex we employ 
the following minmod function = ^xiiU{Xi)\i^j<i): 

=max|min{(t/(a;,)|,+,<i,[/r"},f/r°}, (3.35) 

where we subsequently reset the vertex value to U{xi)\i+j^i := ^^liU {xi)\i+j<i) . 

Proceeding, we estimate the average vertex value over the stencil to its value on the minmod'ed 
neighborhood by computing, Avgg{U{xi)\i+j<i) = | J2e U{xi)\i+j^i, and then we calculate a vertex- 
weighted difference between this average and U^' , which is given by: 

W, = 3 (Avg,(f/(x,)|,+,<i) - f/;^) . (3.36) 
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A schematic of the dynamic adaptive linear restriction method 

• (6>^7^) - vertices of f2e^. 

• U"^^^, ?7°"° - the extrema over the stencil flxi 

• ^ the minmod function at {ie,rii) 

• We - the vertex- weighted difference of averages 

• ^£ - the redistribution factor 

• U{xe)\i^j<i - the updated solution at {^t^rji) 

• Uij\iJ^j<i - the updated monomial coefScents 

(6,^2) 

Base element, Jlg^ £ 



The redistributed slope of 
the solution on 0,^. e Q-^i 




Figure 6: Here we provide a key for the adaptive linear restriction method from §3.4. Again this 
limiting procedure depends on the stencil VLxi) the linear part of the full solution of order p, and 
on a redistribution strategy that can be thought of heuristically as depending on a "consistent 
redistribution of the slopes of the hnear coefficients." 

The restricted difference functions Di are then given with respect to each vertex xi, 

= {U{xe) \i+j<i - U'/) sgnW, (3.37) 

where sgn(-) is the usual signum function except that sgn(O) := 1. Then, if 3)^ is positive, which 
means that either both the average and the approximate solution at the vertex are each larger than 
or similarly that they are both smaller than then we set: 

V — max 1, ^ 1 j , where / = ^^sgnS)^, for each x^ restricted such that 2)^ > 0. (3.38) 
\ m=0 / e 

This allows us now to generate a vertex- wise redistribution factor Me over each element, defined 
simply by setting 

I 0, otherwise, 
where the maximum allowed value is determined by: 

{Uixe)\^+,<^-U^) if sgnW^, > 0, 
{Up'^-U{xe)\i+j<i) otherwise. 



3 A dynamic-in-p family of slope limiters 
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The approximate values at the vertices are then updated, where we make sure the maximum 
redistribution amount is not exceeded, = mm{^i, ^^^^) . The redistributed vertex value is 
updated explicitly to satisfy: 

Uixi)\i+j<i := Uixe)\i+j<i - ^sgnW,. (3.41) 

As an optional step, we add the ability to adapt our limiter to sense areas where substantial 
overshoots and/or undershoots have occurred, thus marking the presence of potential shock fronts. 
We check back to determine that if the redistribution at a specific vertex passes a given tolerance 
e e M"*", then we either zero out the higher order terms if in a fixed order p solution, or we lower our 
polynomial order from p to pum (where pum may be p — 1 or Pmin, e.tc.) if in a p-adaptive context 
(which will be fully addressed in §5). That is, we define a restriction function = [/|j+j>i) 
that operates either on the restricted solution f/|j+j>i or the local polynomial order ^^(figj over 
the entire cell: 




UU+jyi = if iU{x,)\,+,<, - <e)A ^'^^>\Qe,), 

^H^e,) ^ ^''-H^e,) if iUixi)\i+,<, - $a,J < e) A (j9-adaptive) A ^'+^>\n,^) 



if any of the vertex values exceed the tolerance. We also note that clearly e should have an implicit 
dependence on h. 

Finally we make sure that the difference is properly re-weighted for the next computation at 
the elements next vertex (if one exists) by determining the amount available to redistribute by 
computing: We := (Wi — ^esgnWe). This proceeds until no vertices are left to evaluate in the cell. 

Thus we arrive with the sub-quadratic approximate solution, but what we need are the coeffi- 
cients on (f)iQ and 0oi in the basis. To get these we must simply invert the following local constant 
matrix: 

000(3:1) (j)loix2) (l)0liX3)\ /Uoo\ /U{Xi)\i+j<i\ 

(j)oo{xi) (j)io{x2) (poiixs) p^io = U{x2)\i+j<i , (3.42) 

<Poo{Xl) (pw{X2) (pOliXs)/ \Uoi/ \U{x3)\i+j<iJ 

which provides the unknowns. 



§3.5 The hierarchic Unear recombination 

Now, we develop a new slope limiting strategy based on the limiter presented in §3.4, but trans- 
formed into the Taylor basis <jjj G J^(]R^, A^), and generalized over linear recombinations of linear 
reconstructions. 



More clearly, we take our transformed solutions (3.16) such that in the Taylor basis we can 
extract the hierarchical basis at any level [, independently of cell vertices Xi, by simply extracting 
for any hierarchical index b the set {^b,^h+g,'^b+g+i} from §3.2. Notice that this set is entirely 
determined by its indices i and j by way of b{i,j). That is, we can simply denote {^6, '^b+g, ^b+g+i} 
as the first three coefficients of the [i + j)-th derivative of U^. As in §3.2 this provides our linear 



reconstruction, such that equation ( 3.21[ ) becomes our effective sub-quadratic restriction of the 



{i + j)-th derivative of U]^ which we substitute into the formalism of §3.4. That is we set 



U{xi)\i_i>+j_j/<l —^h + '^b+giVi ~ Vc) + '^6+9+1 (6 ~ ^c), 
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A schematic of the hierarchic hnear recombination method 



{^e,^i) - vertices of figj 

b{i,j) - the monomial indices 

^max^ jjrnm _ ^j^^ exticma ovci tlic steucU VtXi 

^ the minmod function at {^e,rii) 
We - the vertex-weighted difference of averages 
^e - the redistribution factor 
U{xe)\i-i'+j-j'<i - the updated 

linear recombination at {C,e,Vi) 
f/(j_i/)(j_jv)|i_j/+j_j/<i - the updated 

monomial coefficents at level i{i + j) 




Base element, Qe £ 



The redistributed slope of 



the solution on Qp G Q 



Figure 7: Here we provide a key for the hierarchic linear recombination method of §3.5. This 
procedure depends on the chosen stencil Qxi, a collection of linear recombinations of restricted 
subsets of monomial coefficients from the total solution, and an application of the method developed 
in 53.4 to these linear recombinations in order to recover the limited solution. 



where i' and j' correspond to the sub-quadratic polynomial basis in the derivation of f/^ with 
coefficients at level i{i + j); or, correspond to the coefficients of the linear recombination at level 

Then (3.34) is calculated, where we evaluate over every in decreasing order. That is, for 
b + g + 1 < s, we compute starting at the top {k — l)-st order derivative steps (3.34)-(3.42) from 
§3.4 with respect to each base coefficient b at that level l{k — 1). Then, due to the redundacy 
of representation for the mixed terms as discussed §3.3, we employ any of our minmod functions 
$J from §3.3 (note that in the experiments below we always use the MUSCL minmod function). 
This is performed until we reach the level corresponding to 6 = 1, at which point we perform the 
calculation one more time identically to that presented in §3.4 except in the Taylor basis. 

Notice here that when the top order is linear, or when p = k = 1 the strategy from §3.4 is 
equivalent to §3.5 up to a change of basis (for example in (3.42) the 
provides for identical error behavior at p = = 1. 



's become ^ij's), which 



§4 Slope limiting: numerical results 

In this section we solve two example problems for an advected scalar quantity l = L{t,x). All of 
our solutions have been run in parallel using an upwinding scheme for the choice of flux. 
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§4.1 Convergence of solutions 

The examples developed in §4.2 and §4.3 both display discontinuities that have meaningful affects 
on the theoretical rates of convergence. Thus first we simply restrict to a smooth solution. That is, 



we use the same formalism of a scalar transport equation (4.1) developed in detail in §4.2, though 



in this case we change the initial conditions to a smooth Gaussian centered at the origin, given by 
lq = age"^^ +^ where = 1 and the boundary condition is the standard transmissive condition 
on both Lf, and Ui,. This is a steady state Gaussian field that "rotates" about the origin by way of 
a pseudo-timestepping. The convergence results are shown in Table [T} Figure |8] and Figure [9j 



/7-convergence 




h level 



Figure 8: The regression rates of convergence for the p G {1, . . . , 5} cases are given by the slope of 
a linear regression line taken from the data in Table [1} 



p-convergence 




Figure 9: The convergence in p G {1, 
Table [H 



5} for the different mesh sizes, as taken from the data in 



§4.2 The rotating half annular crest, cone, and hill solution 

Here we solve a standard rotating landscape solution to a scalar transport equation. That is, 
consider the hyperbolic advection problem: 

dtL + u-V^L = 0, (4.1) 
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p 


LV^°^-crror 


Lf^^ projection error 


X = 1/h 


1 


1.76 X 10-75.31 X 10-^ 


4.57 X 10-^ 


64 


2 


3.34 X 10-74.75 X 10-^° 


1.63 X 10-11 


64 


3 


6.24 X 10-173.19 X 10-11 


8.61 X 10-15 


64 


4 


1.36 X 10-171.63 X 10-12 


8.22 X 10-19 


64 


5 


4.14 X 10-176.54 X 10-1^ 


5.31 X 10-22 


64 


1 


1.05 X 10^76.84 X 10^^ 


7.12 X 10-^ 


32 


2 


3.94 X 10-71.10 X 10^^ 


1.07 X 10-9 


32 


3 


1.55 X 10-73.79 X 10-1° 


2.16 X 10-12 


32 


4 


5.35 X 10-172.49 X 10-11 


8.52 X 10-16 


32 


5 


1.67 X 10-11/1.12 X 10-12 


2.12 X IQ-i^ 


32 


1 


5.24 X 10-76.17 X 10-^ 


1.07 X 10-^ 


16 


2 


4.44 X 10-72.78 X 10-^ 


6.85 X 10-^ 


16 


3 


3.56 X 10-72.31 X 10-^ 


5.12 X 10-1° 


16 


4 


2.55 X 10-75.73 X 10-1° 


8.78 X 10-13 


16 


5 


1.50 X 10-74.36 X 10-11 


7.89 X 10-15 


16 


1 


2.39 X 10-75.17 X 10-^ 


1.30 X 10-^ 


8 


2 


4.73 X 10-76.89 X 10-^ 


3.54 X 10-6 


8 


3 


7.36 X 10-77.58 X 10-^ 


9.55 X 10-^ 


8 


4 


1.21 X 10-76.41 X 10-« 


5.98 X 10-1° 


8 


5 


1.32 X 10-76.07 X 10-9 


2.19 X 10-11 


8 



Table 1: We show the convergence results for the h and p levels whose errors are bounded by 
machine precision after 64 timesteps. The Lf^^ projection error into the basis is also included, 
though, as is clear, these errors are often below machine double precision (~ 1.11 x lO-i^), and 
hence not particularly meaningful. 

with initial-boundary data given by 

t'\t=o = ''O, and tb = 0, 

corresponding to vanishing boundary data, given a time-independent velocity vector field u = u{x) 
with the transported scalar quantity l — L{t,x) in dimension two, such that x — {x,y) and u — 
{u,v). 

We choose a simple square domain Q = [ — |, |]^, with velocity field u = {y, —x). Then letting 



^4 Slope limiting: numerical results 27 



P 



Limiter type 



BJ limiter ' 



53.2 



DEO limiterlE] 
VertexE51l2a.§3-2 



Recombination^^'^ 



BJ limiterEI'§3-2 



Vertex |2lll2a,§3.2 



Recombination 



i3.5 



BJ limiterEl'§3-2 

Vertex E11!2S,§3.2 



Recombination 



i3.5 



BJ limiter^'§3-2 



Vertex mm,^3.2 



Recombination^^'^ 



L error 
L°° error 



1.5x10" 



0.73 



1.1x10"^ 



0.71 
1.1x10-^ 



0.73 



5.9x10" 



0.67 



2.3x10" 



0.74 



2.3x10" 



0.73 



2.2x10" 



0.74 



2.6x10-3 



0.72 

2.6x10-3 



0.72 



2.6x10-3 



0.72 



2.8x10" 



0.72 



2.9x10" 



0.72 



2.9x10"* 



0.72 



Limiter type 



Adapted BJ 



i3.2.1 



BDS limiter[3'§3.4 
Adapted vertex^^-^-^ 



ReconstructiouMuscL n'2a,§3.3 



Restriction H'^^-^ 



Adapted vertex 



j3.2.1 



ReconstructiouENO '^'^'^^■^ 



Restriction Cl'^^-^ 
Adapted vertex^^'^--*^ 



ReconstructiouENO [^'23,§3.3 



Restriction El'^^-^ 



Adapted vertex 



i3.2.1 



ReconstructiouENO [^'23,§3.3 



L error 
L°° error 



1.5x10-3 



0.73 



5.9x10""' 



0.67 
1.0x10-3 



0.72 



5.9x10" 



0.67 



5.0x10" 



0.64 



2.3xl0"3 



0.74 



2.3x10" 



0.73 



4.7x10" 



0.73 
2.5x10-3 



0.73 



2.2x10" 



0.75 



4.8x10" 



0.69 



2.8xl0"3 



0.72 



2.6x10" 



0.73 



Table 2: We give the and L°°-errors of the approximate solutions after one full rotation with 



respect to (4.1 ), setting h = 1/256, At = 1 x 10 ^ and using Runge-Kutta SSP(5, 3). The error ratio 
for the solution with no limiter at p = 1 is L'^/L°° = 2.55 x IO-VO.61, at p = 2 is L'^/L°° = 2.28 x 
10-^/0.44, at p = 3 is LVL°° = 1.71 x 10~V0-35, and for j9 > 3 is unstable. Though, as expected, 
the error in the stable unlimited solutions concentrate along the discontinuities demonstrating sharp 
(> 10% cell-wise in l) overshoots and undershoots. 

To = ^/4 and defining the auxiliary variables 

Ox = X cosTo — y sinTo and Cj^ = y cosr^ + xsinro, 
we take initial data satisfying: 



1, if A, 

to = ■{ 1 — Ba^^, if S < a, 
Wl + cosvrr), otherwise, 



(4.2) 



where 



and 



A = {ao<B<a)A {O^ < ai) , B = \j[Ox--j + O^, 
-'mm (^a, ^ Ol + {Oy + 1/4)2^ , 



r = a 
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Figure 10: Here we show the p = 1 results from Table [2] after one full revolution. The upper left 
is the exact projection at p = 1, the top right is the DEO limiterK^, the middle left is the 
vertex limitert^'2Hl,§3.2^ ^^^q middle right the BDS limiterl^'^^'^, the bottom left the adapted vertex 
limiter^^'^'^, and the bottom right the BJ limiter^l'^^'^. 
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taking a = 0.18, Qq = 0.025 and Oi = —0.23. 

The exact solution may be determined by noticing that since for any F{x,y), where x = x{t) 
and y = y(t), that 



dF „ „ fx' 



which implies that for 



w = ( ^ j = ( ) , we have the system x = y and y' = —x, 



that may be solved by recombining such that the solution to the second order ODE, y" + y = 
can be viewed as a generator of the rotation matrix R about the origin. That is, we obtain the 
clockwise transformation 

R=M -^^"/V (4.3) 
\^smr cose J 

such that Rx yields the exact solution. 

For our numerical experiments, we follow a similar case to that presented in p5j setting our 
mesh width to h = 1/128 and At = 1 x 10~^ in keeping with the CFL condition on hyperbolic 



transport {e.g. see [3H]) ). Let us briefly discuss the results shown in Figure 10 and Figure 11 and 



Table [2] We note that we have run all of our experiments on a regular structured triangular grid. 



In Figure 10 we see the results for linears. The Durlofsky-Engquist-OsherC^] limiter, the vertex 
limiter [25] and the adapted vertex limiter seem to show qualitatively similar behaviors. The Barth- 
Jespersen limiter [6] is slightly more diffuse here at linears (where the adapted Barth-Jespersen 
shows only slight improvement over the native Barth-Jespersen limiter as well), while the BDS 
limiter [7] described in §3.4 shows by far the best L^-error behavior and clearly maintains the best 
signature behavior of the solution everywhere but at the points of discontinuity, where these values 
are tightly redistributed. As previously suggested [25] the vertex limiter and the Barth-Jespersen 
limiter are both quite sensitive to mesh geometries, where the former is better suited in some sense 
to geometries with "sharp angles," and the latter (the Barth-Jespersen limiter) is well-suited for 
regular structured meshes {e.g. Delaunay triangulations) . However, because of the so-called "blind 
diffusion" of both these regimes caused by local extrema — as discussed in §3.2.1 — this behavior 
is not entirely predictable or monotone with respect to mesh regularity, as we see below. The 
hierarchic linear recombination from §3.5, and the hierarchical reconstruction from §3.3 are both 
equivalent by construction to the BDS limiter at p = 1. 

When p > 1 we see an immediate and substantial degradation in the limiting behavior of all 
the regimes, with the single exception of the linear restriction of the BDS limiter from §3.4. This is 
immediately prevalent at p = 2, where the hierarchic linear recombination method from §3.3 is the 
next best limiting regime and yet has an L^-error more than four times that of the linear restriction. 
In fact, the hierarchical reconstruction method from §3.3 may be the most natural extension of the 
BDS limiter to order p, where the choice of linearization is the most direct application of the 
BDS scheme in the Taylor basis. But even here, where at p = 2 we have added only three more 
degrees of freedom to the polynomial hierarchical basis, we see that performing the limiter on the 
linear reconstructions — which amounts to performing the limiting procedure on only two more 
components {i.e. the linear components which are limited with respect to their respective slopes) 
— shows a substantial loss locally in the sharpness of the resolution along the discontinuities. 
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Figure 11: Here we show the p = 2 to p = 4 results from Table |2] after one full revolution. The 
left column shows the linear restriction of the BDS limiter-^ '^^'^ in descending order, while the right 
column shows the next best limiter, in descending order, i.e. at p = 2, p = 3 and p = 4 the 
hierarchical reconstructiouENO 
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The reason for this loss of resolution is not entirely mysterious or unexpected, though previous 
work has demonstrated geometries where this degradation is not immediately observable at 
p = 2, and this behavior seems related to the mass lumping strategy previously discussed in §3.2.1 
(which deserves closer analysis). Nevertheless, here we see that as p increases the number of 
applications of the limiter to the solution increases as a function of the degrees of freedom at the 
(p — l)-st degree {i.e. (p — l)(p — 2)/2). In fact this is true for each of the limiting regimes, with the 
exception of the hierarchical reconstruction methods from §3.3, which actually perform yet another 
iteration of the limiter by employing one of the minmod functions at top order. However, the 
hierarchical reconstruction methods also seem to benefit from the fact that they utilize information 
coming from nonlinearities present in the solution at every level by linearizing with respect to these 
nonlinearities {e.g. equation (3.30)) — in contrast to the vertex-based schemes which linearize about 
a single monomial component (3.21 ) of the expansion, and then utilize a regularizing constraint such 
as (3.24). It turns out that the addition of this nonlinear signature behavior at higher order seems 
to allow the hierarchical reconstruction methods to capture the profile more completely, even with 
the additional application of the limiting regime at each timestep. 

However, by far the most effective limiting regime for p > 1 is the linear restriction of the BDS 
limiter from §3.4, where in 9^ the e has been set to 10~^. Again, this result is not entirely unexpected, 
since slope limiting, as its name suggests, finds its roots in limiting the slopes of lines with respect 
to some linear basis |10]. That having been said, it then seems unlikely that one should be able 
to expect an improvement in the accuracy of a solution near a sharp front simply by applying the 
slope limiter more frequently to the linearization of its respective monomial components. Since, for 
example, if one assumes (fairly realistically) that the top order component has an approximately 
fixed order error which is introduced upon application of the limiter to the FEM solution, then each 
subsequent application of the limiter to the lower level I components should only be able to increase 
the subsequent error introduced over all. In the hierarchical reconstruction methods of [H |27], on 
the other hand, the componentwise minmod function attenuates this effect somewhat, as does the 
fact that all of the limited higher order components serve to help limit the lower order components 
at every level t. 

Before discussing this further, let us first confirm that this result is not simply a special case of 
(4.1) which demonstrates a pathological behavior with respect to (4.2). Below we take a solution 
with admits a number of additional types of singular submanifolds that help to further explicate 
each limiter's behavior. 



]4.3 Steady state convective torque 



Now we show a steady state solution to equation (4.1 ), which effectively isolates the error present in 



the form of torque away from the steady (discontinuous) state in a rotating constant frame solution. 
Our goal here is to present a more difficult set of singular submanifolds Bi <zB present with respect 
to a steady state solution (where the solution here is thought of as the base manifold B) in order to 
more completely isolate the error explicitly introduced by the limiting regimes over varying order 
p. 

Here we work over the Cartesian domain Vt = [0, 2] x [—1, 1], given the same boundary conditions 
from §4.1, and where the exact steady solution is characterized by a velocity field satisfying u = 
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{y, 1 — x) and a steady state scalar field l given by: 

2 



where 



-r 
3' ' 

2ai (1 + cos [(r - 
3ai, 

3a3 (1 + cos [(r - 
ai, 
0, 



a2)vr]) 
a2)vr]) 



1 

4' 



if r < ai 
if ai < r < 3.5a3 
if 4a3 < r < 2ai 

603 < r < 7a3 
if 803 < r < 9a3 
otherwise 



02 



13 



and 



as 



1 

To' 



(4.4) 



P 




error 
L°° error 




error 
L°° error 


1 


BJ limiter^l'^^-^ 

J J fj ±±X±±± UV./X 


f 1.2xl0"2\ 
y 0.49 J 


VprtPvl2Sll2H],§3.2 


( 1.2xl0"2\ 
Y 0.51 J 


1 


DEO limiterES 


f l.OxlO-^A 
0.47 J 


BDS limiterEl'§3-^ 


( 6.8xl0-3\ 
0.40 ) 


1 


Recombination^^'^ 


( 6.8xl0-3\ 
0.40 ) 


ReconstructionENO t^f^'^^-^ 


( 6.8xl0-3\ 
0.40 ) 


2 


BJ limiterEl'§3-2 


( 1.9xl0"2\ 
\^ 0.50 ) 


Restriction El'^^-^ 


( 6.6xl0-3^ 
\^ 0.38 ) 


2 


Vertex mm,^^-^ 


( 1.9xl0-2\ 
\^ 0.50 ) 


Adapted vertex^^-^'-*^ 


( 1.9xl0-2\ 
\^ 0.50 ) 


2 


Recombination^^'^ 


( 1.9xl0-2\ 
Y 0.50 ) 


ReconstructionENO nEZ],§3.3 


( 1.8xl0"2\ 
y 0.52 ) 


3 


BJ limiter-fi-'§3-2 


( 2.2xl0"2^ 
\^ 0.50 ) 


Restriction El'^^ '^ 


( 7.7x10-3 A 
\^ 0.39 ) 


3 


Vertex mm,^^-^ 


( 2.3xl0"2\ 
0.50 ) 


Adapted vertex^^-^--*^ 


( 2.2xl0-2^ 
0.50 j 


3 


Recombination^^'^ 


( 2.3xl0-2\ 
Y 0.50 ) 


ReconstructionENO n'2a,§3.3 


( 2.2x10-2 A 
Y 0.51 ) 


4 


BJ limiter-6 '§3.2 


( 2.3xl0~2\ 
0.50 ) 


Restriction El'S^ "^ 


( 7.7xl0-3\ 
\^ 0.38 ) 


4 


Vertex mm,^^-^ 


( 2.3xl0"2\ 
0.50 ) 


Adapted vertex^^-^'-*^ 


( 2.3xl0"2\ 
0.50 J 


4 


Recombination^^'^ 


( 2.3xl0-2\ 
\^ 0.50 ) 


ReconstructionENO ni'23,§3.3 


( 2.2xl0"2\ 
y 0.51 ) 



Table 3: We give the and L°°-errors of the approximate solutions after T corresponding to a 1/4 
rotation with respect to (4.4), setting h = 1/128, At = 5 x 10"'^ and using Runge-Kutta SSP(5, 3). 
The error ratio for the solution with no limiter at p = 1 is L'^/L°° = 3.9 x 10~3/0.40, at p = 2 is 
L2/L°° = 2.9 X 10-70.34, at p = 3 is LVL°° = 2.4 x IO-VO.23, and for p > 3 is unstable. Again, 
as in Table 2, the unlimited solutions are dominated by local overshoots and undershoots along the 
discontinuities. 



The solution B as shown in Figure 12 is augmented from the relatively well-behaved circular 
convection case analyzed in [25]. Here we have similar outer rings (though substantially "thinned"), 
but have supplemented a pair of inner ring submanifolds that have a thickness of no more than a 
single point that similarly intersects an inner cone along a line of singular points, and with a very 
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thin island outer ring. These initial conditions are not particularly well-behaved, as can be seen 
in Figure 12, where even in the L^-projected exact solution at p = 7 there are variations (jagged 
lines) at the mesh resolution along the lines of singular points. To compound this, we use a larger 
domain than that of |25], which effectively doubles the velocity of the pseudo-timestepping in the 
y-direction, providing for even more instability in the solution space. 

Note that in Figure 12 and Figure [13} the asymmetry in the solution is merely due to that 
fact that we have only gone a quarter turn, thus the diffusive signature of each limiter has only 
been advected a quarter turn, and accumulates or dissipates according to the local behavior of the 
advective flux. 

Now, notice that the adapted limiters from §3.2.1 are not well-suited to handle (4.4) at all. In 
fact (3.26) is, in particular, adapted to represent a case which almost always leads to problems, 
since it does not deal differentially with the special case of U^f^^ = f/,™'°, which in (3.26) up to 
the resolution h is the case for nearly every element in the domain, leading to an almost globally 
uniform "blind diffusion." In fact the adapted cases are almost identical to the native cases at low 
p — when not explicitly dealing with U^^^^ = t/™" — even though the native vertex and Barth- 
Jespersen limiters do not recognize local extrema at all, while the adapted cases do recognize local 
extrema up to, but not including, the degenerate case of t/j™^"^ = ^j™"- As p increases the repeated 
iterations of the limiter swamps this behavior in both the native and adapted limiters, and thus the 
solutions converge to the same value. It is possible that a mass lumping strategy might mitigate 
some of these affects (see [2S] for more information on this technique). 

Moreover, in this example (4.4) the Barth-Jespersen limiter is clearly initially more diffuse than 
the native vertex limiter, which is primarily due here to the fact that the singular submanifolds 
are chosen such that they — again up to the mesh resolution h — spatially oscillate on a local 
neighborhood which is larger than the characteristic length of the edge neighborhood, and so the 
focal neighborhood is a more appropriate area to "sense" in order to capture this semi-localized 
signature behavior. Moreover, the problem of local extrema as discussed in §3.2.1 is of lesser 
importance in this case, since up to the set of codimension one submanifolds of Qh in (4.4), the 
entire domain is characterized and dominated by extremely sharp profiles, making the diffusion — 
which is potentially "blind" near smooth regions — more appropriate here. However, again as p 
increases this behavior gets swamped by the repeated iterations. 

The linear restriction of the BDS limiter ^^'S^ "^ once again demonstrates the best limiting behavior 
as a function of increasing p, which again seems to emphasize the fact that limiting a solution for 
p > 1 must somehow account for the implicit nonlinearity present internal to the cell in a relatively 
explicit way; or, at least, a way which is fully functionally coupled to the entire solution as it exists 
everywhere on the local cell. 

Nevertheless, the linear restriction still substantially outperforms all of the competing limiting 
regimes. There seems to be some indication here that, at least presently, one may expect that near 
areas dominated by shocks the best accuracy that one can hope for is linear accuracy, while still 
hoping to preserve physically important characteristics of the solution {e.g. positivity perserving, 
local conservation of mass, etc.). Of interest, is that this observation falls very neatly in line with 
the state of the art in /ip-adaptive numerical schemes, where a general heuristic follows that for 
potentially discontinuous solutions, in areas of high cell-wise variability, the local order of p is only 
increased if inter-element jumps are small or bounded and controlled, and the internal cell-wise 
variation is strictly bounded above by the cell(s) (usually a subset of cells) containing the global 
maximum [TUl ISO] . 
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Figure 12: Here at top we show the L^-projection of the exact solution at p = 7, with the xz-plane 
slice on the right after 1/4 turn. The middle shows the p = 1 case of the linear restriction El>§3-4^ 
and the bottom shows the p = 1 DEO limiterl^. 
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We explore this issue some in the subsequent section as it applies to p-enrichment, though 
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we also note that at present we are not aware of any formal results which come anywhere near to 
formulating a theorem that subsumes this observational fact (which may in general prove to be only 
one part of the story). Nevertheless, such a result would be of substantial importance to the field, 
as would a counter example, which here could simply be the development of a fully p convergent 
slope limiting regime that limits at all levels I while still preserving the important physical features 
of the solution (and of course does so without relying on prior knowledge, such as the existence of 
an exact solution). 



§5 Adjoining the dynamic /^-enrichment 

Here we present a number of generalizable p-enrichment /de-enrichment schemes based on local 
data and apply them to the problems from §4. These p-enrichment schemes may be viewed as 
alternatives to, for example, the specific energy methods presented in [30] which rely upon the 
variational global entropy of the system of equations, and those discussed in [11] and [5l|19], which, 
as in [30] , try to maximally enrich the domain based on global solution behavior taken with respect 
to the available computational resources and either a priori or a posteriori estimates. 



A general approach based on local data 

We implement a dynamic p-enrichment scheme that utilizes a number of different methodologies in 
order to capture higher order structure in areas of "permissible variability." This scheme is built 
with respect to our collection of p-adaptive slope limiters from §3, such that we inherently arrive 
with a dynamically limited p-enriched solution. 

The nuance of implementing such a scheme in the generalized formulation is that the solution 
must demonstrate a minimal smoothness condition in areas of p-enrichment, while in areas ap- 
proaching discontinuity, p-enrichment must be suppressed in order to maintain stability (especially 
in the absence of a limiter). This issue is not a concern of course when one is able to make smooth- 
ness assumptions a priori about the entire solution space over Q x (0,T) {viz. the formalism of 
[9j and [23j), and has been shown to demonstrate very nice behavior especially in solution spaces 
which are not only smooth, but where in particular one would like to resolve stable areas of maximal 
variation {e.g. as are applicable in some storm surge model applications |23j). 

Nevertheless, in the context of a slightly more generalized system of equations with, for example, 
a coupled hyperbolic equation (or possessing a hyperbolic character in a system of equations) such 



as (4.1), such assumptions cannot generally be made over the entire discrete solution space over 
flh X (0,T), since areas demonstrating strong local gradients VxUh may indicate the presence or 
formation of numeric shock fronts (even given smooth initial data), in which case local p-enrichment 
has a destabilizing effect on the solution (that is, the weak approximation to a discontinuity becomes 
more ill-behaved with respect to increasing p). 

Here we are concerned with dynamically p-adapted solutions to the generalized formulation of 



(2.13) and (2.14) in conjunction with the slope limiters presented in §3. We implement a very simple 
set of p-enrichment strategies, which as we will see, generally tend to undersample the variational 
space {e.g. in contrast to, for example, the poor man's or poor man's greedy algorithm of [TT] 
which always adapts based on some percentage of a global relative bound). The reason for this 
simplification here is to reduce the number of varying parameters in the scheme, in order to isolate 



^5 Adjoining the dynamic p-enrichment 37 



the stabihty of the solution with respect to the hmiting schemes of §3. Hence, we simply set hard 
tolerances which do not depend on, for example, the available computational resources or global 
bounds on the solution. 

Now, in order to additionally deal with both smooth and discontinuous initial-boundary data 
(as well as smooth and discontinuous solutions in (0,T)) we implement the following two distinct 
dynamic p-enrichment schemes — namely we designate them: Type I and Type II p-enrichment 
schemes. We also note that in this section all functions are defined with respect to the master 
element Ai representation. 

The first type of enrichment scheme {i.e. Type I ) applies to solutions in which smoothness may 
be assumed a priori over the entire domain Q x (0,T). That is, taking the approximate solution 
vector Uh we compute the auxiliary sensor over each i-th component of of the state variable U 
(having m components, as in §2): 



U' 



u 



h\c 



Xj 



(5.1) 



where c is the centroid of element f2e and ooj is the midpoint of the j-th edge of fig, and the solution 
Uh is evaluated at these two points, respectively. For smooth solutions, the function Xj may be set 
to either the distance Xj = ~ c| as in [23], or the product Xj = as in [S]. In either case, over 
each timestep n the following p-enrichment functional = ^eX^^i^ei)) evaluated over each 
cell Vt(,{. 







Type I p-enrichment 






if ((sup, sup^. nj > e) A (A; + 1 < p^,J) V (tq > r) , 






if (infi sup^. n;. < e) A (A; - 1 > p^^) A (tq > t^), 






otherwise. 




(5.2) 



where tq is a counter that restricts the p enrichment /de-enrichment such that it may only occur 
every V" timesteps, and where G {1, . . . ,p}. 

For solutions demonstrating approximately nonzero local approximate gradients VxUh 7^ 0, 
wherein we might expect local discontinuities we must find an estimate of the local relative "smooth- 
ness" of Uh- One way of doing this is by setting the auxiliary sensor equal to the following Van 
Leer minmod function across elements (as used in [9]): 



where Vj is the j-th vertex of 



minmod (i/^j I ^ 
As W 



UWcUll: 



Ul\c), (5.3) 
— )■ the solution becomes smoother, and one may 



subsequently employ (5.2) 



A slightly simpler method of dealing with discontinuous solutions simply using local information 
is to define a local smoothness estimator (as discussed in [JT] and [33]) such that we again may 
calculate an elementwise version of (5.1) depending only on the the interior of figp such that: 



\ui-u 



h\\Ll{Qei) 



\u 



h\\Li{nei) 



(5.4) 
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p 


Limiter type 


Type I, L2 


Type II, L2 


e 


c, 


c 




q 


1-5 


BJ limiterlS'§3-2 


3.66 X 10-3 


2.18 X 10-3 


0.1 


-1 


0.1 





2 


1-5 


VertexESl^.§3-2 


3.10 X lO--'^ 


2.10 X 10-3 


0.1 


-1 


0.1 





2 


1-5 


Restriction ^'^^■'^ 


X 


8.03 X 10-^ 


0.1 


-1 


0.1 





2 


1-5 


Recombination^'^'^ 


2.17 X 10-3 


2.09 X 10-3 


0.1 


-1 


0.1 





2 


1-5 


ReconstructionENO '^'^'^^•^ 


1.98 X 10-3 


1.91 X 10-3 


0.1 


-1 


0.1 





2 



Table 4: We give the L^-errors of the approximate solutions after T corresponding to a 1/4 rotation 
with p-enrichment on (4.4), setting h = 1/128, At = 5 x 10~^ and using Runge-Kutta SSP(5,3). 



for the L'^ norms (except when g = 2 in which case we take the standard inner product, as used in 
our examples below), where IJh is the elementwise projected solution ^^~^{Q^), such that in our 



mixed version (5.2) becomes: 





Type II p-enrichment 




if (sup, logio nr < A) A (A; + 1 < p^ax), 




if (infi logio >A)^{k-l> p^in) A (ro > r). 




otherwise, 


(5.5) 



where the bound satisfies 



A 



logio ck 



for p > Pmin 

otherwise 



(5.6) 



SUPi logio ■ 

such that c, c G are user defined constants, where c G (0, 10) is recommended (see for exam- 
ple [4T]) for resolving discontinuities in the context of /ip-adaptivity, and where we have found 
c G (—2,2) optimal. The basic intuition that underpins the use of (5.4) is the observation that 
discontinuous basis functions are assumed to decay, for smooth solutions, at a rate comparable to 
that of the Fourier coefficients in a standard expansion of the solution — which clearly decay at 
a rate of 1/fc^ for q = 2 (see EH HI]); to which we obtain an indicator of the relative local 
regularity of the solution, i.e. the faster the coefficients decay, the smoother the local solution. 



Thus we obtain equation (5.4), which approaches zero as the solution becomes smoother, where 
setting c > is a sharper restriction than the more permissive [i.e. less stable) condition c < 0. 

The results are shown in Table |4] and Figure 14 As expected from before, the linear restriction 
from §3.4 is again by far the most accurate of the choice of limiters when it is stable, where it is 
important to note that in the p-enrichment case the restriction function 9i from §3.4 is calculated 
using e = 10~^, which has the effect of passing cells containing steep gradients to the dynamic 
p-enrichment functions This is enough, it turns out, to make the Type I p-enrichment regime 
unstable with respect to the dynamically adaptive linear restriction limiting regime from §3.3 due in 
part to the function of 5^, which creates an unstable p-fiickering along sharp profile edges. However, 
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Figure 14: Here we show the p values mapped over the p = 1, L^-projected profiles at T = 0.003 
on the top solutions, and at T = 0.06 on the bottom solutions using the settings given in Table 
H The solut ions on the left use the Type II p-enrichment, and those on the right use the Type I 
p-enrichment. 



even turning off the p-de-enriching functionality of 9^ does not help in this case, since the linear 
restriction still zeros out the higher order components, which in the Type I case effectively still 
allows p to flicker locally, leading to the formation of instabilities along sharp edges. We also note 
that in Table |4] we have suppressed the L°°-error, as numerical experimentation suggests that very 
small changes in the p-enrichment settings e, c, c, and can cause big shifts in L^^, which make 
the L°°-error a deceptive measure in the discontinuous p-enrichment case. 

Finally, we emphasize that the p-enriched slope limited solutions show substantially better ac- 
curacy than the constant-in-p solutions from §4.2. This can be attributed in large part to the 
observation that the majority of error in the solutions is accumulated along the discontinuities, 
which is precisely where the p-enrichement schemes p-transition between levels (see Figure 14). 
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Hence, in both Type I and Type II cases {i.e. spatially, from either side of the discontinuity) the 
p-enrichement strongly attenuates (by explicit truncation) the oscillatory instabilities present in 
these regions, as long as the solution does not flicker unstably between them. 

§6 Conclusion 

We have presented a discontinuous Galerkin finite element method for solving dynamically p- 
enriched solutions with consistent slope limiting to arbitrary order in two spatial dimensions over 
generalized coupled systems of PDEs. We have provided a formalism for transforming between 
the polynomial basis of different regimes in order to move between representation spaces. We then 
introduced, up to but not including a substantial choice of minmod functions, seven dynamic-in-p 
slope limiting regimes, and performed numerical experiments on these regimes in order to develop 
a sense of their strengths and weaknesses. We found that our numerical results suggest that, given 
discontinuous initial data, slope limiting over fixed order solutions when p > 1 is most effectively 
accomplished by restricting back to the linear case and using a sharp limiter in that regime, rather 
than keeping the higher order data and trying to limit it in a consistent way — which we found 
introduces more numerical diffusion {i.e. error) on average over time. 

We then presented two types of p-enrichment schemes, fully coupled to the above slope limiting 
regimes. These schemes are designed to exploit certain properties of the solution, and simple 
algorithms were implemented. We then tested these coupled systems on the same model problem 
in order to develop a sense of how dynamic-in-p systems perform relative to fixed-in-p systems. 
Here again, we found that restricting to the linear case seems to be the most effective (and also, 
incidentally, efficient) way of limiting a dynamically p-adapting solution. Moreover, we found that 
in general using the Type I and Type II methods of p-enrichment the accuracy of the solution was 
substantially improved {i.e. by an order of magnitude) with respect to the native solution using 
only the dynamic-in-p slope limiters of §3. 

Future directions include taking the slope limited solution from §3 coupled to the p-enrichment 
scheme from §5 and adding dynamic /i-adaptivity to it, in order to fully exploit the power of hp- 
adaptive convergence. 
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